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MINIMUM-VARIANCE  STATE  ESTIMATION  FOR  UNIFORM  CAUSAL  FUNCTIONAL  EQUATIONS 
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Abstract 


Recently,  a new  class  of  time-domain  state 
space  models  has  been  developed  (Ref.  11  and  12) 
to  describe  layered  media  systems.  When  layers 
are  uniform,  the  resulting  state  equations  are 
referred  to  as  uniform  causal  functional  equa- 
tions. In  this  paper  we  develop  the  minimum- 
variance  state  estimator  for  such  equations.  We 
are  led  to  a natural  form  of  parallel  data  pro- 
cessing for  practical  implementation  of  our 
estimator. 


1.  Introduction 

During  the  past  few  years  there  has  been 
great  Interest  In  modeling  layered  media  systems 
(Refs.  1-12,  for  example)  in  such  diverse  areas 
as  reflection  seismology,  transmission  lines, 
speech  processing,  optical  thin  coatings  and 
EM  problens.  Traditionally,  the  models  have  been 
transfer  function  models;  however,  recently  Nahl 
and  Mendel  (Ref.  11)  and  Mendel,  ct  nl.  (Ref.  12) 
have  shown  that  such  systems  can  be  modeled  In 
the  tine-domain  by  what  appears  to  be  a new  class 
of  state  equations,  causal  functional  equations. 
These  equations  arise  in  a natural  way  when  model- 
ing lossless  layered  media  which  arc  described  by 
the  wave  equation  and  boundary  qondltions,  through 
the  use  of  ray  theory. 

Causal  functional  equations  are.  In  general, 
linear  cent! nuous-time  equations  with  multiple 
tine  delays  (due  to  layers  of  different  travel 
times).  They  do  not  contain  integrals  or  deri- 
vatives; hence,  they  are  not  Integral  or  differen- 
tial equations;  nor  are  they  finite-difference 
equations.  As  is  the  case  with  delay-time  systems, 
they  require  initial  value  information  over 
Initial  intervals  of  time.  Because  of  their  pure 
delay-nature,  their  impulse  response  is  comprised 
of  an  infinite  sequence  of  non-uniformly  spaced 
impulse  functions  ! pen  Mendel,  et  al.  (Ref.  12) 
for  specific  details  into  the  nature  of  causal 
functional  equation). 

In  this  paper  we  direct  our  attention  at  the 
special,  but  very  important  and  useful,  case  of 
uniform  causal  functional  equations,  in  which  all 
tine  delays  are  equal.  This  occurs  when  all 
layers  have  equal  travel  times  (layers  of  unequal 
travel  Line  can  be  built  up  from  layers  of  equal 
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travel  time  by  use  of  Interfaces  with  zero  reflec- 
tion coefficients).  We  shall  develop  a minimum- 
variance  state  estimator  for  uniform  causal  func- 
tional equations. 

Our  uniform  causal  functional  state  space 
model,  denoted,/  in  the  sequel,  is  described  by 
the  following  state  and  observation  equations: 

x(t+x)  - A x(t)  + B m(t)  + w(t)  (1) 


J. 


*U)  - H x(t)  + n(t)  (2) 


mcRr  Is  a known  Input,  wcRn  is  a ran- 
dom disturbance,  £tRs,  ncRs  is  random  measurement 


In  yi  , xcR  _ _ 

ance,  ycRs,  ncRs  is  random  aeasur 
noise,  AeR"*",  BeR**1-,  and  HcR**n.  For  J , the 
following  initial  Interval  Information  is  assumed 
known: 


x(o)  where  o e and  ^ * (0,x) 


o> 


In  Eqs.  (1)  and  (3),  x Is  the  uniform  time  delay. 

An  example  which  Illustrates  the  genesis  of  Eqs. 

(1)  and  (2)  for  seismic  waves  In  layered  media  is 
given  In  Reference  11. 

Some  discussion  Is  in  order  regarding  noise 
processes  w(t)  and  n(t) . For  differential  systems 
ve  usually  assume  that  w(t)  and  n(t)  are  white 
noise  processes.  In  which  case  ETw(t)  w’(()}  ■ 

Q 4(t-{)  and  E(n(t)  n’(C))  - R 4(t-£).  For  such 
systems,  the  correlation  value  must  contribute 
finite  values  in  infinitesimal  intervals.  In  order 
to  have  nonzero  measure,  the  correlation  value  roust 
be  infinite  over  the  infinitesimal  interval.  This 
problem  does  not  arise  for  causal  functional  equa- 
tions. In  fact,  we  would  not  want  v(t),  for  ex- 
ample, to  be  a white  noise  process ^for,  then  g(t) 
would  also  be  white  noise,  which  does  not  make 
sense. 

We  shall  show,  in  Theorem  1 below,  that  causal 
functional  equations  have  solutions  that  are  quite 
similar  to  the  solutions  of  flnltc-differcnce 
equations.  Our  second-order  noise  statistics 
should  therefore  be  closer  to  those  of. a discrete- 
time  system.  (Recall  thar  If  w(k)  and  n(k)  are 
discrete-time  white  noise  sequences,  then  E(w(l) 
w’(J)>  - Q4j j and  E{n(t)  n'(j))  - R4ij  where  4,j 
Is  the  Kronecker  delta  function,  which  Is  equal  to 
unity  for  i-j  and  Is  equal  to  zero  for  i^J.) 

DtSfRBOTION  STATEMENT  A 
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In  out  work  wo  assume  chat  w(t)  and  n(t)  arc 
autually  uncorrelated  Ornstein-Uhlenbeck  processes 
(kef.  16  and  17).  According  to  Parzen  (Ref.  18), 
a stochastic  process  (x(t),  t>0)  is  said  to  be  an 
Ornsteln-Ulilenbcck  process  with  parameters  a>0 
and  B>0  If  it  is  a gausslan  process  satisfying 
E(x(t)l  • 0 and  F.(x(t)  x(C)  l * a e“8|C-t|.  Feller 
(Ref.  19)  states  that  the  Ornstcln-Uhlenbeck 
process  is  obtained  by  subjecting  the  particles 
of  a Brownian  motion  to  an  elastic  force. 

Parzen  states  further  that  the  Ornsteln-Uhlen- 
bcck  process  is  a model  for  Brownian  motion 
which  is  somewhat  more  realistic  than  a Wiener 


We  can  achieve  the  effect  of  a Kronccker 
delta  in  the  second-order  statistics  of  the 
Ornstcln-Uhlenbeck  process  by  letting  B-***.  In 
this  case,  E(x(t)  x(O)  • a If  t - £ and  E(x(t) 
X(()>  ■ 0 if  t 4 £.  More  succinctly,  for  an 
Ornstein-Uhlenbeck  process  (x(t)>  t>0)  for  which 
6-**,  E(x(t)  x(£) ) " a its'  We  shall  assume  that 
w(t)  and  n(t)  are  vector  Ornstein-Uhlenbeck 
processes,  each  of  whose  components  has  the  pre- 
ceding property;  l.e.,  we  assume  that  w(t)  and 
n(t)  are  gausslan  processes,  for  which 

E(w(t) ) - 0 and  E{n(t)>  - 0 V tcR  ft  [0,-)  (4) 

E(w(t)  w’(£))  • Qitc  Vt,  (tH  (5) 

E(n(t)  n'(£)>  ■ Ritr  * t,  £c*  (6) 


E(w(t)  n'(£) ) 


V t,  £cH 


To  complete  the  description  of  J we  must 
specify  the  statistics  of  our  initial  Interval 
Information  in  (3).  We  shall  assume  that  x(o) 
(aij)  is  also  an  Ornstein-Uhlenbeck  process  for 
which  fl— , and  that  x(o)  is  uncorrclated  with  w(t) 
and  n(t);  l.e.,  we  assume  that  x(o)  (oc  j ) Is 
gausslan,  and 


E(x(o) ) 


V oij 


E(x(t)  «’(£))  - A4tC  V t.  £ev? 


E(x(t ) w’(£)>  - o and  E(x(t)  n'(O)  - 0 

V tej  and  tc®.  (10) 

In  many  layered-media  systems  (e.g.,  a lay- 
ered earth  system)  x(o)  • 0,  V ot\J  . For  those 
systems  A * 0 and  the  two  conditions  in  Eq.  (10) 
are  satisfied  because  x(o)  - 0,  V oc\f  . For  the 
sake  of  generality,  we  present  our  results  below 
for  arbitrary  x(o),  oev  . For  > he  sake  of 
simplicity,  all  our  results  are  given  for  time- 
invariant  and  stationary  systems.  It  is  straight- 
forward to  state  them  for  time- varying  and  non- 
statlonary  systems  (l.e.,  for  the  case  when 
A - Aft).  B - B(t) , H - H(t) , Q - Q(t),  R - R(t), 
and  A • A(t)]. 

An  we  have  already  mentioned,  system  A is  a 
conti nuous-time  system  that  is  closely  related  to 


a discrete-time  formulation.  If,  for  example,  we 
set  t - kt  In  4 , whore  k - 0,  1,  2,  ....  and  we 
assume  that  m(t)  and  w(t)  are  only  available  at 

sample  points  t » 0,  r,  2t then  A reduces  to 

an  equivalent  discrete-time  system.  Usually,  how- 
ever, r Is  much  larger  than  real  data  sampling 
rates  (in  geophysical  applications,  data  is  common- 
ly sampled  at  1 or  2 msec  rates);  hence.  It  would 
be  necessary  to  insert  many  small- layers,  whose 
Interfaces  have  zero  reflection  coefficients,  to 
convert  A into  a practical  discrete-time  system. 
That  system  would  be  of  very  large  dimension,  and 
it  Is  doubtful  that  a Kalman  filter  for  this  high- 
order  system  would  be  practical.  For  example.  If 
T - 20  msec  and  dntals  sampled  every  1 msec,  each 
layer  (which  is  described  by  2 states,  an  upgolng 
and  a downgoing  state)  would  be  replaced  by  20 
layers  that  would  be  described  by  2 * 20  - 40 
states.  A 100  layer  system  would  be  described  by 
4,000  states.  By  our  approach,  that  same  100  layer 
syatem  would  be  described  by  200  states. 

Section  2 presents  some  important  preliminary 
results  pertaining  to  the  solution  of  Eq.  (1)  and 
the  statistics  of  x(t)  and  £(t).  A minimum-vnr lance 
state  estimator  is  derived  In  Section  3 first  for 
the  case  when  B « 0 and  then  for  the  more  general 
case  when  B i 0.  Discussions  on  an  Implementation 
for  our  state  estimator  are  given  in  Section  4. 


2.  Preliminary  Results 

In  this  section  we  present  results  pertaining 
to  the  solution  of  Eq.  (1)  and  the  statistics  of 
x(t)  and  jr(t).  Throughout  this  section  and  the 
rest  of  the  paper  we  use  the  fact  that  tc.fi.  can  be 
uniquely  characterized  by  the  mapping 

t • t ' + Mr  where  t'cj  and  M Is  an  Integer  (11) 

We  depict  this  mapping  In  Figure  1. 

Theorem  1.  The  solution  to  the  uniform  causal 
functional  equation* 

x(t+t)  - A x(t)  + w(t)  ; x(o)  otJ  (12) 


xlt*  + (k+l)r)  - Ak+1  x(t’)  + l Ak-1  w(t’+ir) 

1-0 


where  t'c^  , k • 0,  1,  2,  ....  and  t - t'  + (k+l)r. 


Observe  that  (13)  explicitly  shows  how  the 
state  at  any  time  t - t'  + (k+l)r  depends  on  an 
Initial  condition  x(t')  and  the  input  w.  It  is 
of  Interest  to  note  that  x(t)  depends  only  on  a 
single  element  of  the  initial  values  x(o)  (oc«?  ) 


*Our  preliminary  results  are  given  In  this  section 
for  the  B - 0 case.  Their  extensions  to  the 
B i 0 case  are  not  needed  for  our  Section  3 
results. 


:!ion  f) 
ion  □ 

□ 


lOISI.VC 

Dill 


I 


naaely  x(t'),  and  a finite  number  of  point  values 
of  w.  This  shows  that  the  solution  to  the  uniform 
causa  1 fund  1 onalst  ate  equation,  a lthough  con- 
tinuous-! liw-^ln  nature,  derives  Its  values  In  a 
disc r e te-tlnie  fashion  for  a Riven  fixed  value  of 
t'ty^  . Ot  course,  there  are  an  uncountable  num- 
ber of  polncs  In ^ ; hence,  we  can  Imagine  x(t) 
as  being  generat ed  by  an  uncountable  number  of 
discrete-time  systems. 


E(x(t1)x,(t2)} 


E(x(t '+Mt)x* (t"+Jx) J 


E{[AMx(t’) 
(AJx(tn)  + 


M-l  Mil 

+ l A"'1  w(t'+lT)) 
1-0 

J-l  ill 

l AJ_I_J  w(t"+Jt)]M 


j-0 


(20) 


Proof  of  Theorem  1:  Me  use  an  Iterative 
argument  to  justify  the  solution  of  Eq.  (12), 
given  In  Eq.  (13).  First  pick  a t'c*/  . From 
Eq.  (12),  we  see  that 


If  t*  4-  tM,-  then  the  right-hand  side  of  Eq.  (20) 
evaluates  to  zero,  because  of  Eqs.  (5),  (9),  and 
(10).  We  must  now  ascertain  when  t'  i t". 

He  shall  show  that 


x(t’  + t)  - A x(t’)  + w(t')  (14) 

Again,  from  Eq.  (12),  we  sec  that 

x(t’  ♦ 2t)  - A x(t*  + t)  + w(t'  + t)  (15) 
Substitute  Eq.  (14)  Into  Eq.  (15)  to  show  that 

x(t'  + 2t)  - A2x(t')  + Aw(t')  + w(t'  + r)(16) 


t’  i t"  Iff  |t1-t2|  - kt  k - 0.1.2,/..  (21) 

Our  proofs  of  both  the  necessity  and  sufficiency 
of  Eq.  (21)  are  by  the  method  of  contradiction. 

To  begin,  we  assume  the  truth  of  t*  i t"  and  assume 
also  that  - tj  - Pt,  where  P - ... ,-2,-1 ,0,1,2, 
...  . From  Eqs.  (19a)  and  (19b),  this  means  that 

t’  - t"  + (M  - J)t  - Pt  (22) 


Iterating  Eq.  (12)  In  this  same  manner  k+1  times, 
ve  obtain  the  solution  form  in  Eq.  (13). 


Corollary  1.  Both  x(t)  and  £(t),  te(2  , are 
zero  mean  gaussian  process. 

Proof : By  assumption,  x(o)  (oe^)  and  v(t) 
(tefi.  ) are  gaussian  processes.  Equation  (13) 
demonstrates  chat  x(t)  (t c<8>)  can  be  decomposed 
Into  a linear  combination  of  these  processes; 
hence,  since  linear  combinations  of  gaussian  pro- 
cesses are  also  gaussian  (Ref.  13),  x(t)  is  gaus- 
sian V tcC.  Since  x(o)  (ac\f  ) and  w(t)  (te(ft  ) 
are  by  assumption  zero  mean,  it  follows,  from 
Eq.  (13)  that  x(t)  (tc  (2  ) is  also  zero  mean. 

Because  the  proof  of  the  results  for  £(t)  is 
so  similar  to  that  Just  given  for  x(t),  we  leave 
It  to  the  reader. 

Next,  we  present  some  important  cross-covari- 
ance results. 

Theorem  2.  Let  t.  and  t.  be  points  In 

Then 


Efxftj)  x’(t2))  - 0 

(17) 

,|  * kT  k - 0,1,2,... 

(18) 

This  theorem  states  that  x(t)  is  uncorrelated 
at  nonlntcger  multiples  of  the  uniform  delay  time 
t. 

Proof : Time  points  tj  and  tj  can  be  express- 
ed by  means  of  Eq.  (11)  as 

tj  - t'  ♦ Ht  t'tV  and  M an  Integer  (19a) 

tj  ■ t"  ♦ Jt  tMew?  and  J an  Integer  (19b) 

From  Eqs.  (19a),  (19b)  and  (13),  we  find  that 


t*-t"“(P-M  + J)t  (23) 

Now  t’,  t "cm/  - (0,t);  hence, 

-T  < t'  - t”  < T (24) 

Consequently,  since  P,  M,  and  J are  all  Integers, 
the  only  way  Eq.  (23)  can  be  true  is  if  P-M+J  « 0, 
in  which  case  t'  - t";  but,  this  contradicts  the 
assumption  that  t'  i t".  He  conclude,  therefore, 
that  If  t'  j*  t",  then  tj  - t2  d Pt,  P • ...,-2,-1, 

0,1,2 which  can  also  be  stated  as  Uj-tjl  i 

kT,  k - 0,1,2,...  . 

Next,  we  assume  the  truth  of  U^-tjl  I * kt, 
k - 0,1,2,...,  and  assume  also  that  t'  - t".  From 
Eqs.  (19a)  and  (19b),  this  means  that 

tx  - t2  - (M  - J)t  ; (25) 

but,  this  contradicts  our  assumption  that  1 1 j^— t j i »* 
kt.  We  conclude,  therefore,  that  if  Itj-tjl  t ki, 
then  t'  i t".  This  completes  the  proof  of  Thoerem 
2. 

Corollary  2.  Let  t^  and  t2  be  points  In  . 

Then 

ECiUjjx.'dj))  - 0 (26) 

and 

E<*<t2)*,(tl)}  " 0 • (27) 

If 

|tj  - t2|  i kT  k - 0,1,2,...  (28) 

Proof:  Froa  Eq.  (2),  we  see  that 


Since  t'  i t",  we  conclude,  f roe  Eq.  (21) 


I{X(t1)i,(t2>)  - HEUa^x'CtjjjH* 

+ E(n(tj)n'(t2)  ) + } 

♦ B(£(tI)s,(t2))N'  (29) 

The  flret  tern  In  this  expression  is  zero  via 
Theorem  2.  Observe  that  Eq.  (28)  implies  that 
*1  ^ l2<  hence,  the  last  three  terms  In  Eq.  (29) 
arc  xero  via  Eqs.  (16),  (13)  and  (10).  The  proof 
of  Eq.  (27)  follows  in  a similar  manner. 

In  Section  3,  where  we  derive  the  mlnlmum- 
varlance  estimator  for  x(t),  we  assume  that  mea- 
surements, /(t),  are  available  from  time  zero 
up  to  and  Including  time  t.  We  now  introduce 
three  data  sets.  Let  Y(t)  be  the  set  of  mea- 
surement data  which  are  available  up  to  and  includ- 
ing tine  t;  l.e., 

T(t)  - (y(l):  tc£>  (30) 

Let  Y|f(t')  be  a finite  set  of  measurements  associa- 
ted with  point  t 'cs?  and  Integral  multiples  of  x 
from  t';  l.e., 

TM(t’)  - (x(f+jT),  J - 0,1,2 M)  (31) 

For  J - 0,  Yfj(t')  only  contains  the  single  measure- 
ment x^'j*  E°r  j • H,  Y(((t *)  contains  Mel 
measurements  made  at  the  points  labeled  t'.  t'+t, 

. ...t'+Mt  in  Figure  1.  Finally,  let  Yc(t')  be  the 
set  difference  between  Y(t)  And  YM(t');  l.e., 

Yc(f)  - Y(t)  - YM(t*)  (32) 

Observe  from  Eqs.  (30),  (31),  and  (32)  that 


T(t)  - YM(t')  U Yc(t') 

(33) 

and 

Vfjnyf)  - ♦ 

(34) 

where  p Is  the  null  set. 

Theorem  3.  Let  t be  a point  In  & such  that 
the  representation  In  F.q.  (11)  Is  valid.  Then 
Y|.(t')  Is  statistically  Independent  of  Yjft')  and 
x?t)  Is  statistically  Independent  of  Yc(t'). 

Proof:  We  represent  t as  In  Eq.  (11) 
choose  a tvpical  element,  X^l2)»  fro™  data 
Yc(f);  l.e.. 

and 

sec 

X(t2)  t Yc(t') 

(35) 

By  construction  of  Yc(t')  we  know  that 

t2  - t"  ♦ Jr  j » 0,1,. ...M 

(36) 

and  t"c,J  . We  observe  that  t"  i t';  for, 
t'  * t"  then  t2  • t*  + jr  (J  * 0,1,...,M)S 

If 

and 

these  time  points  are  associated  with  YM(t'). 
From  our  preceding  discussions,  we  know  that 
YM(t')  and  Yr(t')  share  no  comon  time  points; 
hence,  t”  i t'. 


that 

|'t  - t2|  i kr  k » 0,1,...  (37) 

which  means  that  time  points  In  Yc(t')  satisfy  Eq. 
(28);  hence,  from  Eq.  (27),  we  see  that 

E{x(t2)x’(t))  - 0 for  -»i(t2)  e-Vc(t*)  (38) 

Since  x(t)  and  x<0  are  gausslan  (Corollary  1), 
the  fact  in  Eq.  (38),  that  x(t?)  and  are  un~ 

correlated  for  V X^2^  E ¥c(t')  *1®°  means  that 
X(t2)  and  x(t)  are  statistically  independent  (Ref. 
13)  for  V x<t2>  e Yc(t'). 

Next,  pick  a typical  element,  x(l,)»  f rom 
YM(t’);  l.e., 

X(tx)  c YM(t')  (39) 

By  construction  of  Yu(t'),  we  know  that 

t2  - t'  + jr  J - 0,1,2 M (40) 

Since  t'  i t",  we  know  that  t,  and  t2  satisfy  Eq. 
(28);  hence.  It  follows  from  Eq.  (26)  that 

E{x(t1)x’(t2)}  ■ 0 V x(t2)  c Yc(t’) 

and  X(t2)  c Y^t')  (41) 

From  Corollary  1,  this  once  again  means  that  x^lj) 
and  x(e2)  are  statistically  Independent  for 
V xCt?  e vc(f)  and  x(ti)  t Y^(t ' ) ; or,  that 
Yc(t')  and  YM(t')  are  statistically  Independent. 
This  completes  the  proof  of  Theorem  3. 


3.  Minimum-Variance  State  Estimator 

To  begin,  we  present  the  minimum-variance 
atate  estimator  for  system  /l  with  B»0.  The  case 
of  a known  forcing  function  in  the  state  equation 
la  treated  below  in  Theorem  5. 

Theorem  4.  For  system  J , with  B=0,  the 
minimum-variance  estimator  of  x(t),  where  t * t' 

+ Mr,  t'e J and  M an  Integer,  is 

x(t)  - E(x(t)|YM(t’))  (42) 

This  theorem  states  the  interesting  result 
that  the  optimal  estimator  of  x(t)  need  not  be 
conditioned  on  the  entire  data  set,  Y(t).  The 
optimal  estimator  of  x(t)  need  only  be  conditioned 
on  data  set  Y^(t'),  a data  set  which  contains  only 
a finite  number  of  points* 

Proof;  From  estimation  theory  (Refs.  14  and 
IS,  for  example).  It  Is  well  known  that  the  minima 
variance  estimator  of  x(t)  Is 

x(t)  ■ E(x(t)  [ Y ( t ) ) (43) 

which  can  also  be  written,  using  Eq.  (33),  as 

x(t)  - E(x(t)|YM(t’),  Yc(t’>)  (44) 
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Next,  we  use  the  fact  (proven  In  Medltch  (Ref. 

14) • PS*  101)  that  If  a,  b,  and  c are  gaussian 
random  vectors  and  b and  £ are  ata» latlcally 
Independent,  then 

E(a|b,c)  * E(a|b)  + E{a |jc ) - E{a)  (45) 

Since  jc( t ) , YM( t ' ) and  Yc(t’)  are  gaussian  and 
Yx(t')  and  Y^(t')  arc  statistically  Independent 
(Theorem  3),  ve  can  expand  the  right-hand  side 
of  Eq.  (44)  by  means  of  Eq.  (45),  to  obtain 

x(t)  - E(x(t)|YM(t'))  + E(x(t) |Yc(t')) 

- E{x(t) ) (46) 

Additionally,  since  x(t)  and  Yc(t')  are  statisti- 
cally independent  (Theorem  3),  E(x(t)  |Y<.(t') ) ■ 
E(x(t)};  thus,  Eq.  (46)  reduces  to  Eq.  (42). 

Corollary  3.  For  any  fixed  t't^  , x(t)  * 
x(t'+Mt)  can  be  computed  by  the  discrete-time 
Kalman  filtering  equations  with  t'  considered  as 
the  Initial  starting  time.  Hence,  for  t'cj  , 
compute  x(t'+Mi)  from  the  following  equations: 

xft’+Nt)  - Ax[t ' + (M-l)t)  + K(t'+MT)(£(t’+Mt) 

- HAx[ t ' ♦ (M-l)t) ) (47) 

Plt'+Milt’  + (M-I)tJ  - 

AP(t ' + (M-1)t | t ' + (M-1)t)A*  + Q (48)* 

K(t'+Mx)  - Plt’+«T|t'+(M-l)T) 

• H’(HF[t  *+Mt |t'+(M-l)T)H’  + R)-1  (49)* 

and 

P(t'+Mt|t'-HlT)  - ( I-K(t '+Mt)H) 

• Plt'+HT|t,+(M-l)T)  (50) 

where  M ■ 1,2,...,  and,  x(t')  ■ E(x(t'))  and 
P(t'lt')  - EWt'Jx'd'lT  - A. 

Proof:  Equation  (42)  states  that  x(t)  Is 
obtained  from  a finite  point  measurement  data  set. 
It  can  be  expressed  by  means  of  Eqs.  (11)  and  (31) 

•s 

x(t)  - E{x(t,+Mr)|z(t').jr(t'+t) *(t'+Hi))  (51) 

The  right-hand  side  of  Eq.  (51)  Is  analogous  to 
the  optimal  filtered  estimate  for  a discrete-time 
system  with  sampling  rate  equal  to  t sec  and  Ini- 
tial time  equal  to  t'.  Our  notation  In  Eq.  (47) 
for  x(t'+Mt)  is  equivalent,  therefore,  to  the 
notat  ion  x(t  '+>lT  1 1 ’+Mt)  . Equations  (47)  through 
(AO)  are  the  discrete-time  Kalman  filter  equations 
(Ref.  14,  for  example)  obtained  by  replacing  dis- 
crete variable  k in  those  equations  by  t'  + Mr . 

•For  nonstationary  noise  process,  Q becomes  Q 1 1 * 

♦ (M— 1 ) T ) and  R becomes  R(t'  ♦ Mr). 


He  observe,  from  Corollary  3,  that  for  each 
fixed  t'e^,  the  optimal  estimate  of  the  state 
at  t'  + Mr  (M  - 1,2,...)  Is  obtained  by  Iterating 
the  familiar  Kalman  filter  equations  on  T,  using  t' 
as  the  Initial  starting  time.  Since  there  are  an 
uncountable  number  of  points  In ^ , this  results 
In  an  uncountable  number  of  estimators  which  have 
to  be  Implemented  In  order  to  obtain  x(t)  for  any 
tc  fi.  . .We  discuss  a practical  Implementation  of 
theae  results  In  Section  IV.C. 

Next,  ve  shall  extend  our  Theorem  4 and  Col- 
lary  3 results  to  the  case  where  there  is  a known 
forcing  fucntlon  In  the  state  equation,  that  Is, 
to  the  case  when  BpO  in  . 

Theorem  5.  Let  x(t)  denote  the  minimum-vari- 
ance estimator  of  x(t)  for  /i  In  Eqs.  (1)  and  (2). 
Let  Xj(t)  denote  the  minimum-variance  estimator  of 
x(t)  under  the  assumption  that  B-0  In  Eq.  (1). 
Additionally,  xj(t)  Is  associated  with  the  follow- 
ing modified  measurement  equation: 

Xx(t)  • X(t)  - H i(t)  (52) 

where  |(t)tRn  satisfies  the  uniform  causal  func- 
tional equation 

*(t+r)  » A *(t)  + B m(t)  (53) 

for  which 

l(o)  - 0 V oc  J . (54) 

Then, 

x(t)  « XjU)  + lit)  . (55) 

The  proof  of  this  theorem  la  given  In  Appendix 
A.  The  essence  of  the  theorem  Is  that  x,(t)  can  be 
computed  using  Theorem  4 or  Corollary  3 for  the 
modified  measurement  y^(t)  which  can  be  constructed 
a priori  from  £(t)  and  lit).  Signal  £(t)  can  be 
precomputed  since  n(t)  Is  known  a priori.  The 
solution  of  Eqs.  (53)  snd  (54)  Is  given  by  Eq.  (13) 
where  x(t')  - 0 and  w is  replaced  by  Bm.  He  then 
obtain  the  desired  estimate  of  ait),  x(t),  by  the 
simple  linear  transformation  In  Eq.  (55). 

Theorem  5 provides  the  formal  results  for 
x(t).  We  do  not  recommend  calculating  x(t)  by 
Theorem  5;  for  the  calculations  Involve  solving 
an  auxiliary  causal  functional  equation.  A more 
practical  way  for  computing  x(t)  Is  given  In  the 
following: 

Corollary  4.  Let  x(t)  denote  the  minimum- 
variance  estimator  of  x?t)  for  & In  Eqs.  (1)  and 
(2).  For  any  fixed  t'c J , x(t)  ■ i(t'Hit)  can 
be  computed  from  the  following: 

x(f+MT)  • Ax[t'+(M-1)T)  ♦ Bmlt’+fN-Dr) 

+ K(t'+MT){£(t'+m)  - HAx(t '+(M-1)t| 

- HBm[t '+(M-l)t] ) , (W 

where  F(t '+Mt  1 1 '+(M-l)t J , K(t'+flT),  and 
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P(t '+Mt | t '+Mt)  are  given  by  Eqs.  (48),  (49),  and 
(SO),  respectively,  and  H - 1,2....  . Addition- 
ally, i(t')  - E(x(t ' ) } and  P(t'lt')  - E(x(t ') 
x'(t')T  -A. 

Proof:  For  t * t'+Hr,  we  can  express  x(t'  + 
Nr)  by  aeana  of  Eqs.  (SS),  (S3)  and  (47),  as 

*(t’+m)  - ^(t'+Hr)  + 4(t'+MT) 

- Ax1It,*(M-l)Tj  ♦ KU'+MtH^U'+Ht) 

- HAx^ [t ’+(M-1) r ] ) + A ilt'+fM-Drl 

♦ B5(t,+(M-1)T)  (57) 

where,  froa  Eqs.  (S2)  and  (S3), 
y1(t'+MT)  • i(t'+Mr)  - H j(t'+Mx) 

- jt(t'-rtlT)  - HAgJ t '+(M-l)x ) 

- HBm[ t '+(M-l)x ) (58) 

Substitute  Eq.  (58)  into  Eq.  (57),  making  use  of 

Eq.  (SS),  to  show  that. 

i(t’-Hh)  - Ax( t ,+(M-1)t ] + Bo(t'+(M-l)T] 

♦ K(t  ’♦Mr) (jf(t  '+Ht)  - HAx( t '+(M-1)t] 

- HBm[t '♦(M-l)r ) ) (59) 

which  Is  our  desired  result,  Eq.  (56) . Observe 
that  gain  matrix  K(t’+Mt)  In  Eq.  (59)  is  the  same 
gain  matrix  as  In  Eq.  (47);  hence.  It  Is  compuced 
from  Eq.  (49),  and  P| t '♦Mi j t '+(M-l)t ) and 
P(t  '♦Mr  1 1 ’+HO  are  computed  froa  Eqs.  (48)  and 
(SO),  respectively. 

We  see  that  a known  forcing  function  in  a 
uniform  crusal  functional  state  equation  Is  handled 
In  the  estimator  equation  In  exactly  the  same 
manner  that  such  a function  Is  handled  for  a 
finite-difference  equation.  It  does  not  affect 
the  error  covariance  equations. 


4.  Computation 

A.  Introduction 


B,  Digital  Computer  Simulation  of  Uniform  Causal 
Functional  Equation 

We  assume  here  that  t Is  an  Integral  multiple 
of  sample  rate  T;  l.e.,  we  assume  that 

T * LT  (60) 

The  approximation  denotes  a temporaL. quantization 
which  may  be  needed  to  associate  a value  of  L with 
t.  For  r « LT,^  ♦ $ where 

Jc  - (t * : t'  - IT,  1 - 0,1, . . . ,L-1 ) (61) 

Set  J (■  has  a countable  number  of  elements.  It 
does  not  Include  t-T,  since  set^  does  not  Include 
that  point. 

We  direct  our  attention  at  Eq.  (13)  which  is 
the  solution  to  uniform  causal  functional  equation 
(21).  Observe  from  Eqs.  (60)  and  (61)  that 


t'  + pr  - (t  + I.p)T  . • 

(62) 

Letting 

Kjip)  4 y{(t  + Lp)T] 

(63) 

and 

yt(p)  & w[(l  + Lp)T]  , 

(64) 

We  see  that  (13)  can  be  written  as 

xt(k+l)  - Ak+1  xt(0)  + l Ak_i  w^l) 

(65) 

1-0 


where  1 - 0,1, ...  ,1.-1 . This  equation  Is  the  well- 
known  solution  to  the  following  finite-difference 
equation  (Ref.  20,  for  example): 

x^(k+l)  - A Xj(k)  + yt(k)  (66) 

where  l - 0,1,...,L-1  and  k - 0,1,2, ... ,k*.  We 
conclude,  therefore,  that  the  sampled  solution  to 
our  uniform  causal  functional  equation  can  be 
generated  as  depleted  In  Figure  2. 

To  begin,  we  must  fix  our  useful  data  length 
at  an  Integer  multiple  of  t,  say  Nt . In  terms  of 
our  sampling  rate,  T,  our  useful  data  length  Is 
NET  sec. 


In  Section  2,  we  demonstrated  that  the  solu- 
tion, x(t),  of  causal  functional  equation  (1) 
can  be  generated  by  an  uncountable  number  of 
discrete-time  system.  In  Section  3 we  demonstra- 
ted that  x(t)  (tc  R ) can  be  generated  by  means  of 
an  uncountable  number  of  discrete-t ime  Kalman 
filters.  When  we  simulate  our  results  on  a digi- 
tal computer,  computations  are  made  every  T sec, 
at  discrete  time  points.  Consequently,  on  a 
digital  computer,  x(t)'  is  generated  by  a count- 

Iable  number  of  discrete-time  systems,  and  x(t)  Is 
generated  by  a countable  number  of  dlscrete-t lme 
* Kalman  filters.  We  explore  the  detailed  ramifica- 

tions of  these  observations  below. 


On  a discrete  time-scale  the  Ornsteln-Uhlen- 
beck  process,  with  8-—,  reduces  to  a dlscrctc-time 
white  noise  sequence  (w(JT),  j - 0,1,..., J*}  where 
E{w(JT)  w'(lT)}  - QAj, . In  order  to  cover  the 
useful  data  length,  we  choose 

J*  - NL  (67) 

we  create  sequences  w^fk),  Wj^k),  ....  wL_1(k), 
which  are  needed  to  simulate  Eq.  (66)  by  sorting 
and  distributing  (w(JT),  J - 0,1,..., NL)  according 
to  the  distribution  algorithm  given  on  Figure  2. 

Index  k ranges  from  zero  to  k*.  To  determine 
k*,  we  observe  that  the  very  last  value  of  w[(t  + 
Lk)T]  used  by  the  distribution  algorithm 
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(•et  t - L-l  and  k - k*)  is  y{ [ (k*+l)L-l]T) . 

Since Is  an  open  Interval  on  its  right-hand  end, 
the  argument  of  this  function  must  equal  NL-1; 
hence, 

k*  - N-l  (68) 

Our  L noise  sequences,  Wjj(k),  w^(k),  ..., 
wj_i(k),  drive  L systems,  each  given  bv  Eq.  (66), 
which  operate  in  parallel.  Solving  our  uniform 
causal  function al  equation  on  a digital  computer 
has  led  to  a parallel  data  processing  algorithm. 

As  a final  step  we  generate  (x(jT),  J ■ 0, 

1 NL-1)  by  multiplexing  Xfl(k),  XjM.  ..., 

>i -1 (k)  using  the  multiplexor  algorithm  also 
given  on  Figure  2.  It  Is  straightforward  to 
show  that 

(x(JT),  j-0,1 NL-1)  - (x^O.x^O) x^CO), 

XqU)  .x^l) Xj^d) ....  ,Xg(k*)  ,x^(k*)  , 

....^(k*))  (69) 

Also,  given  k and  1,  j ■ t+Lk. 

C.  Digital  Computer  Simulation  of  Minimum- 
Variance  Estimator 

We  now  direct  out  attention  at  the  simulation 
of  Eqs.  (47) - ( 50) . As  in  Section  B,  we  assume 
that  time  is  discretized  so  that  Eqs.  (60)  and 


(61)  are  true.  Letting 

*t(M|«)  k *1(1  + LM)T]  (70) 

K,(M)  £ K[(l  + LM)T]  (71) 

Pt(M|n-l)  A P[(l  + LM)T|(1  + LM  - L)T)  (72) 
P1(M|N)  4 P((t  + LM)T|  (1  + LM)T]  (73) 

and 

Ft(M)  A y((l  + LM)T]  , (76) 


We  refer  to  Eqs.  (’5)  - (80)  as  the  1th  Kalman 
filter,  KF1. 

We  obtain  x(jT|jT)  (J-1,2 J*  - NL-1)  by 

means  of  the  parallel  data  processing  algorithm 
depicted  In  Figure  3. 

The  Kalman  filters  do  not  generate  estimates 
over  the  initial  time  Interval  -4 c;  for,_over  that 
interval  .Xj[(0|0)  is  given  a priori.  The  first 
measurement  used  by  KFq  Is  jr(LT) , whereas  the  last 
measurement  used  by  KFj_^  1®  X.[(NL-1)T);  hence,  we 
must  have  the  sequence  (£(JT) , j“L,L+l, . . . ,NL-1 ) 
available  at  the  start  of  the  simulation.  Kalman 
filter  outputs  Xfl(M|M) ,x, (M|M) , . . . ,Xl_j(M|M)  are 
multiplexed  to  give  the  desired  estimates 

{x(JT|jT),j-L,L+l NL-1)  - (XjjdlD.Xjdll), 

. . . .x^d  1 1)  ,Xq(2  1 2)  ,^(2 1 2) x^d  1 2)  , 

. . . .x^N-llN-l)  .x^N-llN-l) , 

...,xL_1(N-l|N-l))  (81) 

For  a time-varying  system  J , we  must  make 
the  following  substitutions  in  Eqs.  < 75) — ( 78) : 

A£(M)  -*■  A,  H£(M)  -►  H,  Qt(M-l)  - Q,  and  R£(M)  - R. 
where,  for  example, 

At(M)  A A[(t  + LM)T)  (82) 

In  this  case,  each  of  the  L Kalman  filters  has  no 
common  calculations,  and  the  computational  burden 
can  be  quite  heavy. 

If,  on  the  other  hand,  system  Jl  is  time-in- 
variant or  slowly  time-varying  (i.e.,  all  matrices 
are  piecewise  constant  over  t sec  Intervals),  then 
Eqs.  (76),  (77),  and  (78)  are  not  functions  of  1 
(remember  that  1 is  used  to  define  ^ q,  and  r?  c 
is  less  than  T units  in  length).  Hence,  in  these 
cases  we  need  only  calculate  P(M|M-1),  K(M),  and 
P(M|M)  for  M * 1,2,..., N-l  once.  These  calcula- 
tions are  then  used  by  all  L Kalman  filters,  which 
greatly  reduces  the  computational  burden. 


« 


? 


we  rewrite  Eqs.  (4 7 ) - ( 50)  as: 
xt(M|M)  - AXj(M-l  | M— 1 ) + Kl(M)(Jrt(M) 

- HAxt(M-l[M-l) J (75) 
Pt(M|N-l)  - AP£(M-1|M-1)A'  + q (76) 
K,(M)  - P1(m|m-1)H'(HP1(M|m-1)H'  + R]_1  (77) 


A flowchart  for  implementing  the  L Kalman 
filters,  when  is  time-invariant  or  slowly  time- 
varying,  is  depicted  in  Figure  4.  Because  the 
outer  loop  varies _ M and  the  Inner  loop  varies  1, 
outputs  P(M|M)  and  Xg(M|K)  are  generated  in  multi- 
plexed ordering.  This  can  be  verified  by  listing 
the  sequence  of  outputs  from  this  flowchart  and 
comparing  them  with  the  right-hand  side  of  Eq . (81) 
to  see  that  they  are  identical. 


and 

P{(m!«1  - (I  - Kj(M)H)  PjfMjM-l)  (78) 

where  for  each  C ( (*0, 1 , . . . ,L-1) , M * 1,2,..., 

M*  - N-l , and, 

xt(0|0)  - E(xt(0))  (79) 

Pt(0|0)  - A . (80) 


5.  Conclusions 

In  this  paper  we  have  derived  the  minimum- 
variance  state  estimator  for  uniform  causal  func- 
tional equations.  These  equations  are  useful  for 
modeling  layered  media  systems  whieh  are  described 
by  the  lossless  wave  equation  and  boundary  condi- 
tions. Causal  functional  equations,  though  con- 
tinuous-time In  nature,  bear  a strong  resemblance 
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to  diacret-time  equations.  In  fact,  we  have 
shown  that  the  solution,  x(t),  of  uniform  causal 
functional  equation  (1)  can  be  generated  by  an  un- 
countable number  of  discrete-time  systems.  Ue 
have  also  shown  thac  for  any  fixed  t'e » (0,i), 
x(t),  where  t “ t'  + Mr  (M  « 1,2,...),  is  given 
by  the  usual  discrete-time  Kalman  filter  equations 
with  t'  considered  the  initial  starting,  time.  Of 
course,  to  obtain  x(t)  for  all  tc  (P-  we  would  need 
an  uncountable  number  of  discrete-time  Kalman 
filters;  but,  imposing  a mesh  on>?  leads  to  a 
countable  number  of  Kalman  filters  which  operate 
in  parallel,  as  depicted  in  Figure  3.  To  the 
best  knowledge  of  the  authors,  this  is  the  first 
estimation  theory  result  that  has  led  to  a natural 
form  of  parallel  data  processing. 

The  results  of  this  paper  are  not  merely  an 
end  unto  themselves.  An  important  problem  for 
layered  media  systems  is  to  extract  reflection 
coefficients  from  noisy  measurements.  This  if 
often  referred  to  as  an  inverse  problem  (Refs.  3, 
4,  7,  8 and  10),  and  usually,  solutions  are  given 
only  for  noise-free  measurements.  The  reflection 
coefficients  appear  in  matrix  A.  Ey  means  of  the 
results  of  this  paper,  two  appraoches  can  be 
studied  for  solving  the  inverse  problem.  In  the 
first  approach,  we  augment  scate  equations  (one 
for  each  unknown  parameter)  to  Eq.  (1)  and  develop 
an  extended  minimum-variance  estimator  for  the 
resulting  augmented  system.  Of  course,  the 
augmented  state  equations  must  also  be  causal 
functional  in  nature  or  else,  if  they  were  dif- 
ferential equations,  our  augmented  system  would 
be  hcriditary  in  nature  and  a different  course 
of  action  would  have  to  be  taken.  In  the  second 
approach,  we  estimate  the  reflection  coefficients 
by  a maxi mum- likelihood  technique.  To  accomplish 
this,  we  must  first  develop  the  correct  likelihood 
function  for  a uniform  causal  functional  equation. 
Both  of  these  approaches  are  presently  under 
investigation. 
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Appendix  A.  Proof  of  Theorem  5 

Let  xj(r)  be  associated  with  system  , 
when  B-0;  that  is  to  say,  Xj(t)  satisfies 

XjU+t)  * AXj(t)  + w(t)  (A-l) 


where 

x^ (o)  4 x(o)  V oe-^  (A-2) 

Add  Eqa . (A-l)  and  (S3),  and  (A-2)  and  (S4),  to 
show  that 

Xjlt+r)  + g(t+x)  » AlXjU)  + _g(t)l  + Ba(t) 

+ w(t)  (A-3) 

and 

XjOj)  + g(o)  - x(o)  V oe</  (A-4) 

Comparing  Eqa.  (A-3)  and  (1),  and  (A-4)  and  "(3), 
we  conclude  by  uniqueness,  that 

x(t)  - Xj(t)  + g(t)  (A-5) 

Next,  we  define  the  modified  measurement  vec- 
tor  £,(t)  in  Eq.  (52).  That  ^(t)  can  be  express- 
ed solely  in  terms  of  x^(t)  is  apparent,  when  Eqs. 
(2)  and  (A-5)  are  substituted  into  Eq.  (52);  i.e., 

ij(t)  - HXj(t)  + n(t)  (A~6) 

Observe  that  Eqs.  (A-l)  and  (A-6)  are  now  of 
the  right  form  to  use  Theorem  4 to  obtain  xx(t)  « 
E{xi(t)  | Y,  (t)  ),  where  Yj(t)  - (jr.(X):  0 <_  X <_  t, 
te~A  ). 

From  estimation  theory,  we  know  that 

*(t)  - E{x(t) | Y ( t ) ) (A-7) 

Applying  (A-7)  to  (A-5),  we  find  that 

x(t)  - EtxjlOlYU)}  + t(t)  (A-8) 

since  £(t)  is  deterministic.  We  must  now  prove 
that 

8(^(0  |Y(t))  - E{x1(t)|Yx(t))  - x(t)  (A-9) 

in  which  case  Eq.  (A-8)  reduces  to  the  desired 
result  in  Eq.  (55). 

We  shall  examine  the  sigma  fieWs  generated 
by  y(t)  and  y^(t).  Let  y and  yj  denote 
the  sigma  fields  associated  with  £(t)  and  2j(t), 
respectively;  i.e.. 


A 

- {w:  £(t,w)  ^ a)  acR® 

(A-10) 

and 

K 

" (id:  2.x(t,u>)  a)  acR® 

(A-ll) 

From  Eq. 

(52),  we  see  that  /Jy  can 

be  expressed 

in  terms 

V 

of  ix(t) , as 

{u>:  ^j(t.w)  + Hg(t)  a)' 

acR® ; 

(A-12) 

hence , 

V 

(w:  2.j(t,w)  < » • H£(t)} 

acR® 

(A-13) 

A 
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If  g(t)  In  finite  for  any  fixed  t,  a - Hg(t)  also 
ranges  through  Ks . This  condition  on  g(t)  Is 
satisfied  as  long  as  Bm(t)  Is  finite  for  any  fixed 
t.  Consequently, 


JB  - (si:  jTjft.w)  < Oj ) OjtR3 
which  proves  that  y " -^y  * 


(A-14) 


Since  ^(t)  and  j^(t)  generate  the  same  sigma 
fields,  conditioning  with  respect  to  jr(t)  [or 
Y(t)]  Is  equivalent  to  conditioning  with  respect 
to  2j(t)  [or  Yj(t)J;  hence,  Eq.  (A-9)  is  true, 
and,  as  we  pointed  out  above,  Eq.  (A-8)  reduces 
to  the  desired  result  in  Eq.  (55). 
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